16S Phyloseq

Load required packages.

require(phyloseq)

Loading required package: phyloseq

 require(parallelDist)

Loading required package: parallelDist
Warning: package ‘parallelDist’ was built under R version 4.1.2

require(ggplot2)

Loading required package: ggplot2

  require(randomForest)

Loading required package: randomForest
Warning: package ‘randomForest’ was built under R version 4.1.2
randomForest 4.7-1
Type rfNews() to see new features/changes/bug fixes.

Attaching package: ‘randomForest’

The following object is masked from ‘package:ggplot2’:

    margin

require(Boruta)

Loading required package: Boruta

  require(phyloseq)
require(stats)
  require(dplyr)

Loading required package: dplyr
Warning: package ‘dplyr’ was built under R version 4.1.2

Attaching package: ‘dplyr’

The following object is masked from ‘package:randomForest’:

    combine

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

require(readr)

Loading required package: readr
Warning: package ‘readr’ was built under R version 4.1.2

  require(pheatmap)

Loading required package: pheatmap

require(tidyverse)

Loading required package: tidyverse
Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ──────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ tibble  3.1.6     ✓ stringr 1.4.0
✓ tidyr   1.2.0     ✓ forcats 0.5.1
✓ purrr   0.3.4     
Warning: package ‘tidyr’ was built under R version 4.1.2
── Conflicts ─────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::combine()       masks randomForest::combine()
x dplyr::filter()        masks stats::filter()
x dplyr::lag()           masks stats::lag()
x randomForest::margin() masks ggplot2::margin()

  require(micrUBIfuns)

Loading required package: micrUBIfuns
Loading required package: data.table
data.table 1.14.2 using 1 threads (see ?getDTthreads).  Latest news: r-datatable.com
**********
This installation of data.table has not detected OpenMP support. It should still work but in single-threaded mode.
This is a Mac. Please read https://mac.r-project.org/openmp/. Please engage with Apple and ask them for support. Check r-datatable.com for updates, and our Mac instructions here: https://github.com/Rdatatable/data.table/wiki/Installation. After several years of many reports of installation problems on Mac, it's time to gingerly point out that there have been no similar problems on Windows or Linux.
**********

Attaching package: ‘data.table’

The following object is masked from ‘package:purrr’:

    transpose

The following objects are masked from ‘package:dplyr’:

    between, first, last

Loading required package: ggsci
Loading required package: ComplexHeatmap
Loading required package: grid
========================================
ComplexHeatmap version 2.8.0
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference

If you use it in published research, please cite:
Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
  genomic data. Bioinformatics 2016.

The new InteractiveComplexHeatmap package can directly export static 
complex heatmaps into an interactive Shiny app with zero effort. Have a try!

This message can be suppressed by:
  suppressPackageStartupMessages(library(ComplexHeatmap))
========================================
! pheatmap() has been masked by ComplexHeatmap::pheatmap(). Most of the arguments
   in the original pheatmap() are identically supported in the new function. You 
   can still use the original function by explicitly calling pheatmap::pheatmap().


Attaching package: ‘ComplexHeatmap’

The following object is masked from ‘package:pheatmap’:

    pheatmap

Loading required package: vegan
Warning: package ‘vegan’ was built under R version 4.1.2
Loading required package: permute
Warning: package ‘permute’ was built under R version 4.1.2
Loading required package: lattice
This is vegan 2.6-2

require(circlize)

Loading required package: circlize
Warning: package ‘circlize’ was built under R version 4.1.2
========================================
circlize version 0.4.14
CRAN page: https://cran.r-project.org/package=circlize
Github page: https://github.com/jokergoo/circlize
Documentation: https://jokergoo.github.io/circlize_book/book/

If you use it in published research, please cite:
Gu, Z. circlize implements and enhances circular visualization
  in R. Bioinformatics 2014.

This message can be suppressed by:
  suppressPackageStartupMessages(library(circlize))
========================================

  require(phylosmith)

Loading required package: phylosmith

require(phangorn)

Loading required package: phangorn
Loading required package: ape
Warning: package ‘ape’ was built under R version 4.1.2

Attaching package: ‘ape’

The following object is masked from ‘package:circlize’:

    degree


Attaching package: ‘phangorn’

The following objects are masked from ‘package:vegan’:

    diversity, treedist

Load necessary files from Dada2 processing script.

load("/Users/gordoncuster/Desktop/Git_Projects/BeeBread/Data/16S/DADA2outputs/dada2_outputs_16S_BeeBread_wtree.RData")

Read in metadata

metadata<-read.csv("/Users/gordoncuster/Desktop/Git_Projects/BeeBread/Data/BB_sample-metadata.txt", sep = "\t")
rownames(metadata)<-metadata$X
metadata$Region_Site<-paste( metadata$Region, metadata$Site, sep = "")
metadata$Treatment_Site<-paste( metadata$Treatment, metadata$Site, sep = "")

Convert files to correct format for phyloseq object and then rename samples in OTU table to match the format of the metadata file. Create phyloseq object and root phylogenetic tree for usage with distance metrics (e.g., Unifrac).

#assign object type to merge into phyloseq object
md<-sample_data(metadata)
otu<-otu_table(seqtab.nochim, taxa_are_rows = F)
tax_tab<-tax_table(taxa)
#check sample names of both otu table and metadata
#remove .fastq.gz from otu table names
sample_names(otu)<-str_split(sample_names(otu), pattern = ".fastq.gz", simplify = T) [,1]
#since they match, you can create the phyloseq object
bb16S_orig<-phyloseq(md, otu, tax_tab, fitGTR$tree)

#root tree for distance metrics like unifrac
set.seed(11)
phy_tree(bb16S_orig)<-root(phy_tree(bb16S_orig), sample(taxa_names(bb16S_orig), 1), resolve.root = TRUE)
is.rooted(phy_tree(bb16S_orig))

[1] TRUE

Pre-processing. Rename OTU IDS to something more manageable. We save the ASCV sequences in case we need them later (e.g., BLASTN searches). Create Bacterial and Archaeal phyloseq objects, but first we remove Chloroplasts.

bb16S_arch

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 9 taxa and 80 samples ]
sample_data() Sample Data:       [ 80 samples by 10 sample variables ]
tax_table()   Taxonomy Table:    [ 9 taxa by 7 taxonomic ranks ]
phy_tree()    Phylogenetic Tree: [ 9 tips and 8 internal nodes ]

Remove control samples as there is only as single replicate of each.

ps_wo_control<-subset_samples(bb16S_bac, Treatment != "none")

Examine rarefaction curves. Plateau in most samples under the 5k read mark.

rarecurve(otu_table(ps_wo_control), step=50, cex=0.5)

First, check sample depths for rarefaction. Minimum sample depth is a bit over 21k reads, so we will rarefy at 20k.

set.seed(11)
sort(sample_sums(ps_wo_control))
bb16S_bac_rarefy<-rarefy_even_depth(ps_wo_control, sample.size = 7000, rngseed = 14, trimOTUs = T)
#360 otus were removed during rarefaction. 

#To move forward, we can hellinger transform and go from there.
bb16S_bac_hellinger<-transform_sample_counts(ps_wo_control, function(x) sqrt(x / sum(x)))

#maximum liklihood point estimates
bb16S_bac_ML<- transform_sample_counts(ps_wo_control, function(x) x / sum(x))

#Exploratory analysis of Alpha diversity

Plots of Alpha Diversity by region, treatment, timepoint, and the combinations.

plot_richness(bb16S_bac_rarefy, x = "Region", measures = c("Shannon", "Observed", "Chao1"), color = "Region") + geom_boxplot() + ggtitle("Region")

plot_richness(bb16S_bac_rarefy, x = "Treatment", measures = c("Shannon", "Observed", "Chao1"), color = "Treatment") + geom_boxplot() + ggtitle("Treatment")

plot_richness(bb16S_bac_rarefy, x = "Timepoint", measures = c("Shannon", "Observed", "Chao1"), color = "Timepoint") + geom_boxplot() + ggtitle("Timepoint")

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Shannon", "Observed", "Chao1"), color = "Site") + geom_boxplot() + ggtitle("Site")

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Shannon", "Observed", "Chao1"), color = "Timepoint") + geom_boxplot() + ggtitle("Site by timepoint")

plot_richness(bb16S_bac_rarefy, x = "Region_Treatment", measures = c("Shannon", "Observed", "Chao1"), color = "Region_Treatment") + geom_boxplot() + ggtitle("Region_Treatment")

plot_richness(bb16S_bac_rarefy, x = "Region_Site", measures = c("Shannon", "Observed", "Chao1"), color = "Region_Site") + geom_boxplot() + ggtitle("Region_site")

plot_richness(bb16S_bac_rarefy, x = "Region_Timepoint", measures = c("Shannon", "Observed", "Chao1"), color = "Region_Timepoint") + geom_boxplot() + ggtitle("Region_Timepoint")

plot_richness(bb16S_bac_rarefy, x = "Treatment_Timepoint", measures = c("Shannon", "Observed", "Chao1"), color = "Treatment_Timepoint") + geom_boxplot() + ggtitle("Treatment_Timepoint")

Based upon the site and region graphs, it appears that site contains more differences than region. There are no obvious differences among treatments. The site by timepoint graph shows what appears to be an interaction of before and after treatment and site. Some sites display lower diversity following treatment and others higher, with consistent trends across all sites though.

Statistical testing of alpha diversity metrics. Using ANOVA framework and checking for assumption of normality and correlation in predictors. region and site appear to be associated, so I will drop Region for this exploratory analysis. We can revisit if we decide Region is a more accurate predictor, but it appears that from the graphs above site provides more resolution than region. https://www.pluralsight.com/guides/testing-for-relationships-between-categorical-variables-using-the-chi-square-test

#create alpha diversity table and prep data to include grouping columns
richness_BB16S_bac<-estimate_richness(bb16S_bac_rarefy, measures = c("Shannon", "Observed", "Chao1"))
#add metadata
richness_BB16S_bac$Region<-sample_data(bb16S_bac_rarefy)$Region
richness_BB16S_bac$Site<-sample_data(bb16S_bac_rarefy)$Site
richness_BB16S_bac$Timepoint<-sample_data(bb16S_bac_rarefy)$Timepoint
richness_BB16S_bac$Treatment<-sample_data(bb16S_bac_rarefy)$Treatment
#add combo columns
richness_BB16S_bac$Region_Treatment<-sample_data(bb16S_bac_rarefy)$Region_Treatment
richness_BB16S_bac$Region_Site<-sample_data(bb16S_bac_rarefy)$Region_Site
richness_BB16S_bac$Region_Timepoint<-sample_data(bb16S_bac_rarefy)$Region_Timepoint
richness_BB16S_bac$Treatment_Timepoint<-sample_data(bb16S_bac_rarefy)$Treatment_Timepoint

#test of association between 
chisq.test(richness_BB16S_bac$Region, richness_BB16S_bac$Site, simulate.p.value = TRUE)

    Pearson's Chi-squared test with simulated p-value (based on 2000 replicates)

data:  richness_BB16S_bac$Region and richness_BB16S_bac$Site
X-squared = 216, df = NA, p-value = 0.0004998

Tests of Shannon diversity - significant differences across sampling sites and interaction term of timepoint and site.

#ANOVA and Tukey's PW comparisons
#Run a full 3 way ANOVA with interactions
mod<-aov(Shannon ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)
#residuals look fine. 
shapiro.test(mod$residuals)

    Shapiro-Wilk normality test

data:  mod$residuals
W = 0.98676, p-value = 0.6541

summary(mod)

                         Df Sum Sq Mean Sq F value   Pr(>F)    
Site                      5  7.215  1.4430   9.262 9.77e-06 ***
Treatment                 2  0.314  0.1572   1.009 0.374654    
Timepoint                 1  0.913  0.9125   5.858 0.020684 *  
Site:Treatment           10  2.196  0.2196   1.410 0.215468    
Site:Timepoint            5  4.597  0.9194   5.901 0.000441 ***
Treatment:Timepoint       2  0.104  0.0519   0.333 0.718673    
Site:Treatment:Timepoint 10  2.244  0.2244   1.440 0.202737    
Residuals                36  5.608  0.1558                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Summary of model shows site to be only lower order significant predictor.
#interaction term of site x time point as well. - this was picked up in the exploratory graphs. 

#What about pairwise differences in site
TukeyHSD(mod, "Site")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Shannon ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)

$Site
              diff         lwr          upr     p adj
Ber-Bel  0.6310367  0.14624536  1.115827989 0.0047899
Dal-Bel  0.1515572 -0.33323409  0.636348544 0.9331063
Sla-Bel  0.3992869 -0.08550437  0.884078264 0.1577373
Spr-Bel -0.3178867 -0.80267800  0.166904632 0.3771633
Sum-Bel  0.4602726 -0.02451876  0.945063865 0.0709627
Dal-Ber -0.4794794 -0.96427076  0.005311869 0.0540091
Sla-Ber -0.2317497 -0.71654104  0.253041590 0.7040609
Spr-Ber -0.9489234 -1.43371467 -0.464132043 0.0000138
Sum-Ber -0.1707641 -0.65555544  0.314027191 0.8938504
Sla-Dal  0.2477297 -0.23706159  0.732521035 0.6432093
Spr-Dal -0.4694439 -0.95423523  0.015347403 0.0623627
Sum-Dal  0.3087153 -0.17607599  0.793506636 0.4094025
Spr-Sla -0.7171736 -1.20196495 -0.232382317 0.0010436
Sum-Sla  0.0609856 -0.42380571  0.545776916 0.9989011
Sum-Spr  0.7781592  0.29336792  1.262950548 0.0003416

plot(TukeyHSD(mod, "Site"), las=1, cex.axis = 0.4)

#Several significant pairwise differences to accompany the significant global test for interaction term. We can
TukeyHSD(mod, "Site:Timepoint")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Shannon ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)

$`Site:Timepoint`
                    diff         lwr         upr     p adj
Ber:A-Bel:A  1.098429161  0.30304953  1.89380880 0.0013633
Dal:A-Bel:A -0.209138823 -1.00451846  0.58624081 0.9984586
Sla:A-Bel:A  0.526283471 -0.26909616  1.32166311 0.4885062
Spr:A-Bel:A -0.296741464 -1.09212110  0.49863817 0.9735081
Sum:A-Bel:A  0.318575516 -0.47680412  1.11395515 0.9564626
Bel:B-Bel:A -0.187446459 -0.98282609  0.60793318 0.9994309
Ber:B-Bel:A -0.023802271 -0.81918190  0.77157736 1.0000000
Dal:B-Bel:A  0.324806822 -0.47057281  1.12018646 0.9504208
Sla:B-Bel:A  0.084843969 -0.71053566  0.88022360 0.9999998
Spr:B-Bel:A -0.526478360 -1.32185799  0.26890127 0.4879560
Sum:B-Bel:A  0.414523126 -0.38085651  1.20990276 0.7972938
Dal:A-Ber:A -1.307567984 -2.10294762 -0.51218835 0.0000891
Sla:A-Ber:A -0.572145690 -1.36752532  0.22323394 0.3655966
Spr:A-Ber:A -1.395170625 -2.19055026 -0.59979099 0.0000279
Sum:A-Ber:A -0.779853645 -1.57523328  0.01552599 0.0588502
Bel:B-Ber:A -1.285875620 -2.08125525 -0.49049599 0.0001186
Ber:B-Ber:A -1.122231432 -1.91761107 -0.32685180 0.0010048
Dal:B-Ber:A -0.773622339 -1.56900197  0.02175730 0.0627758
Sla:B-Ber:A -1.013585192 -1.80896483 -0.21820556 0.0039696
Spr:B-Ber:A -1.624907521 -2.42028716 -0.82952789 0.0000013
Sum:B-Ber:A -0.683906035 -1.47928567  0.11147360 0.1498550
Sla:A-Dal:A  0.735422294 -0.05995734  1.53080193 0.0922320
Spr:A-Dal:A -0.087602641 -0.88298227  0.70777699 0.9999997
Sum:A-Dal:A  0.527714339 -0.26766529  1.32309397 0.4844706
Bel:B-Dal:A  0.021692364 -0.77368727  0.81707200 1.0000000
Ber:B-Dal:A  0.185336552 -0.61004308  0.98071619 0.9994878
Dal:B-Dal:A  0.533945645 -0.26143399  1.32932528 0.4670132
Sla:B-Dal:A  0.293982792 -0.50139684  1.08936243 0.9752414
Spr:B-Dal:A -0.317339537 -1.11271917  0.47804010 0.9575958
Sum:B-Dal:A  0.623661949 -0.17171768  1.41904158 0.2498009
Spr:A-Sla:A -0.823024935 -1.61840457 -0.02764530 0.0371469
Sum:A-Sla:A -0.207707955 -1.00308759  0.58767168 0.9985501
Bel:B-Sla:A -0.713729929 -1.50910956  0.08164970 0.1137108
Ber:B-Sla:A -0.550085741 -1.34546538  0.24529389 0.4228520
Dal:B-Sla:A -0.201476649 -0.99685628  0.59390299 0.9988974
Sla:B-Sla:A -0.441439502 -1.23681914  0.35394013 0.7293475
Spr:B-Sla:A -1.052761831 -1.84814147 -0.25738220 0.0024331
Sum:B-Sla:A -0.111760345 -0.90713998  0.68361929 0.9999966
Sum:A-Spr:A  0.615316980 -0.18006265  1.41069661 0.2666870
Bel:B-Spr:A  0.109295005 -0.68608463  0.90467464 0.9999973
Ber:B-Spr:A  0.272939193 -0.52244044  1.06831883 0.9858021
Dal:B-Spr:A  0.621548286 -0.17383135  1.41692792 0.2540072
Sla:B-Spr:A  0.381585433 -0.41379420  1.17696507 0.8679152
Spr:B-Spr:A -0.229736897 -1.02511653  0.56564274 0.9965020
Sum:B-Spr:A  0.711264590 -0.08411504  1.50664422 0.1163969
Bel:B-Sum:A -0.506021975 -1.30140161  0.28935766 0.5464643
Ber:B-Sum:A -0.342377787 -1.13775742  0.45300185 0.9302956
Dal:B-Sum:A  0.006231306 -0.78914833  0.80161094 1.0000000
Sla:B-Sum:A -0.233731547 -1.02911118  0.56164809 0.9959504
Spr:B-Sum:A -0.845053877 -1.64043351 -0.04967424 0.0291416
Sum:B-Sum:A  0.095947610 -0.69943202  0.89132724 0.9999993
Ber:B-Bel:B  0.163644188 -0.63173545  0.95902382 0.9998429
Dal:B-Bel:B  0.512253281 -0.28312635  1.30763291 0.5285090
Sla:B-Bel:B  0.272290428 -0.52308921  1.06767006 0.9860600
Spr:B-Bel:B -0.339031902 -1.13441154  0.45634773 0.9344903
Sum:B-Bel:B  0.601969585 -0.19341005  1.39734922 0.2952384
Dal:B-Ber:B  0.348609093 -0.44677054  1.14398873 0.9220167
Sla:B-Ber:B  0.108646240 -0.68673339  0.90402587 0.9999975
Spr:B-Ber:B -0.502676090 -1.29805572  0.29270354 0.5561351
Sum:B-Ber:B  0.438325397 -0.35705424  1.23370503 0.7375907
Sla:B-Dal:B -0.239962853 -1.03534249  0.55541678 0.9949490
Spr:B-Dal:B -0.851285183 -1.64666482 -0.05590555 0.0271836
Sum:B-Dal:B  0.089716304 -0.70566333  0.88509594 0.9999997
Spr:B-Sla:B -0.611322329 -1.40670196  0.18405730 0.2750341
Sum:B-Sla:B  0.329679157 -0.46570048  1.12505879 0.9453054
Sum:B-Spr:B  0.941001487  0.14562185  1.73638112 0.0096104

plot(TukeyHSD(mod, "Site:Timepoint"), las=1, cex.axis = 0.4)

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Shannon"), color = "Timepoint") + geom_boxplot() + ggtitle("Site x Timepoint")

Tests of Chao 1 - emphasizes importance of rare taxa. Site is significant as are the interactions with site.

#ANOVA and Tukey's PW comparisons
#Run a full 3 way ANOVA with interactions
mod<-aov(Chao1 ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)
#residuals look fine. 
shapiro.test(mod$residuals)

    Shapiro-Wilk normality test

data:  mod$residuals
W = 0.97941, p-value = 0.2856

summary(mod)

                         Df Sum Sq Mean Sq F value   Pr(>F)    
Site                      5  18455    3691  12.137 6.26e-07 ***
Treatment                 2    926     463   1.523 0.231844    
Timepoint                 1   1068    1068   3.512 0.069058 .  
Site:Treatment           10   7481     748   2.460 0.023369 *  
Site:Timepoint            5   8307    1661   5.463 0.000769 ***
Treatment:Timepoint       2   1854     927   3.048 0.059876 .  
Site:Treatment:Timepoint 10   4913     491   1.616 0.141568    
Residuals                36  10947     304                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Site and the interactions with site are significant. 

TukeyHSD(mod, "Site")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Chao1 ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)

$Site
                diff        lwr        upr     p adj
Ber-Bel  14.53134366  -6.887271  35.949958 0.3402017
Dal-Bel -29.38462857 -50.803243  -7.966014 0.0026464
Sla-Bel -29.40859141 -50.827206  -7.989977 0.0026212
Spr-Bel -20.15949629 -41.578111   1.259118 0.0750034
Sum-Bel -18.02599761 -39.444612   3.392617 0.1416952
Dal-Ber -43.91597222 -65.334587 -22.497358 0.0000059
Sla-Ber -43.93993506 -65.358550 -22.521320 0.0000058
Spr-Ber -34.69083995 -56.109455 -13.272225 0.0002999
Sum-Ber -32.55734127 -53.975956 -11.138727 0.0007291
Sla-Dal  -0.02396284 -21.442578  21.394652 1.0000000
Spr-Dal   9.22513228 -12.193482  30.643747 0.7853117
Sum-Dal  11.35863095 -10.059984  32.777246 0.6066951
Spr-Sla   9.24909512 -12.169520  30.667710 0.7835044
Sum-Sla  11.38259380 -10.036021  32.801208 0.6045698
Sum-Spr   2.13349868 -19.285116  23.552113 0.9996463

plot(TukeyHSD(mod, "Site"), las=1, cex.axis = 0.4)

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Chao1"), color = "Timepoint") + geom_boxplot() + ggtitle("Site x Timepoint")

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Chao1"), color = "Treatment") + geom_boxplot() + ggtitle("Site x Treatment")

Tests of Observed - unique taxa in a sample. Site is significant as are the interactions with site.

#ANOVA and Tukey's PW comparisons
#Run a full 3 way ANOVA with interactions
mod<-aov(Observed ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)
#residuals look fine. 
shapiro.test(mod$residuals)

    Shapiro-Wilk normality test

data:  mod$residuals
W = 0.98392, p-value = 0.489

summary(mod)

                         Df Sum Sq Mean Sq F value   Pr(>F)    
Site                      5  14992  2998.4  12.255 5.64e-07 ***
Treatment                 2    574   287.1   1.173  0.32091    
Timepoint                 1   1267  1266.7   5.177  0.02894 *  
Site:Treatment           10   5400   540.0   2.207  0.04025 *  
Site:Timepoint            5   7668  1533.6   6.268  0.00028 ***
Treatment:Timepoint       2   1550   775.1   3.168  0.05404 .  
Site:Treatment:Timepoint 10   3247   324.7   1.327  0.25368    
Residuals                36   8808   244.7                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Summary of model shows site to be a significant predictor as well as the site interaction terms.
TukeyHSD(mod, "Site")

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = Observed ~ Site * Treatment * Timepoint, data = richness_BB16S_bac)

$Site
              diff        lwr         upr     p adj
Ber-Bel  13.500000  -5.712000  32.7120000 0.3031171
Dal-Bel -27.166667 -46.378667  -7.9546667 0.0018429
Sla-Bel -25.250000 -44.462000  -6.0380000 0.0043102
Spr-Bel -18.666667 -37.878667   0.5453333 0.0609710
Sum-Bel -14.250000 -33.462000   4.9620000 0.2488740
Dal-Ber -40.666667 -59.878667 -21.4546667 0.0000032
Sla-Ber -38.750000 -57.962000 -19.5380000 0.0000080
Spr-Ber -32.166667 -51.378667 -12.9546667 0.0001832
Sum-Ber -27.750000 -46.962000  -8.5380000 0.0014162
Sla-Dal   1.916667 -17.295333  21.1286666 0.9996436
Spr-Dal   8.500000 -10.712000  27.7120000 0.7660676
Sum-Dal  12.916667  -6.295333  32.1286666 0.3499569
Spr-Sla   6.583333 -12.628667  25.7953333 0.9043412
Sum-Sla  11.000000  -8.212000  30.2120000 0.5264786
Sum-Spr   4.416667 -14.795333  23.6286666 0.9817483

plot(TukeyHSD(mod, "Site"), las=1, cex.axis = 0.4)

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Observed"), color = "Timepoint") + geom_boxplot() + ggtitle("Site")

plot_richness(bb16S_bac_rarefy, x = "Site", measures = c("Observed"), color = "Treatment") + geom_boxplot() + ggtitle("Site")

Summary of alpha diversity: Rarefying at 20k reads sufficiently captures the diversity of each sample, as demonstrated by the plateau in our rarefaction curves. Site is the most useful predictor for determining differences in alpha diversity. Several interactions with site are present for all metrics of alpha diversity. No clear trends in alpha diversity differences. Some sites display shifts responding to treatment, while others do not. Some sites show increased alpha diversity at timepoint 2 and others the opposite.

#Multivaraite exploratory analysis

Info on strata argument. Briefly, this allows us to restrict permutations to within a group and control for an effect in an analogous way to a random effect.

https://stats.stackexchange.com/questions/350462/can-you-perform-a-permanova-analysis-on-nested-data/350504#350504?newreg=8e389dc237c54d87ab240486b8b0fe33 https://stats.stackexchange.com/questions/188519/adonis-in-vegan-order-of-variables-or-use-of-strata?noredirect=1&lq=1 https://stats.stackexchange.com/questions/459407/permanova-outputs-with-or-without-random-factor https://ichthyology.usm.edu/courses/multivariate/feb_7.pdf

Both with Bray-curtis and Unifrac distances.

#this model explores the effect of Site, treatment and time point and restricts permutations sites. 
adonis2(dist.uf ~  Site * Treatment * Timepoint, data= sd_adonis, strata = sd_adonis$Site, permutations = 1000)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  strata 
Permutation: free
Number of permutations: 1000

adonis2(formula = dist.uf ~ Site * Treatment * Timepoint, data = sd_adonis, permutations = 1000, strata = sd_adonis$Site)
                         Df SumOfSqs      R2       F  Pr(>F)  
Site                      5  0.18272 0.43042 13.9670 0.03896 *
Treatment                 2  0.00278 0.00655  0.5311 0.70629  
Timepoint                 1  0.00551 0.01297  2.1043 0.14885  
Site:Treatment           10  0.06701 0.15784  2.5610 0.01798 *
Site:Timepoint            5  0.03039 0.07159  2.3230 0.05495 .
Treatment:Timepoint       2  0.00543 0.01278  1.0368 0.36763  
Site:Treatment:Timepoint 10  0.03650 0.08597  1.3949 0.20579  
Residual                 36  0.09419 0.22188                  
Total                    71  0.42453 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Based on results, site is the main determinant of microbiome composition (r2 of 0.51 and r2 of 0.38, BC and weighted unifrac, respectively). Time point is also significant but explains much less of the variation (r2 of 0.063 and r2 of 0.021, BC and weighted unifrac, respectively). Treatment is only significant for BC (r2 of 0.024). Both remain significant for BC dissimilarity when site and timepoint are controlled for using the strata argument. However, when site is controled for using weighted unifrac dissimilarity, the only significant predictor is the site x treatment interaction. Next, I will split by time point and examine the effect of each treatment while controlling for site effects using the strata argument (strata = site).

Take away: Site is the strongest predictor of community dissimilarity. This is true for both Bray-Curtis and weighted unifrac distances.

Split by time point as these are before and after treatment.

#and pairwise
pairwise.adonis2(dist ~ Treatment, strata = "Site",  data= sd_adonis, nperm = 1000)

$parent_call
[1] "dist ~ Treatment , strata = Site , permutations 1000"

$CF_vs_CON
          Df SumOfSqs      R2      F  Pr(>F)  
Treatment  1   0.1655 0.02752 0.6227 0.01299 *
Residual  22   5.8483 0.97248                 
Total     23   6.0139 1.00000                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$CF_vs_ORG
          Df SumOfSqs      R2      F Pr(>F)
Treatment  1   0.1229 0.01954 0.4383 0.1818
Residual  22   6.1706 0.98046              
Total     23   6.2935 1.00000              

$CON_vs_ORG
          Df SumOfSqs      R2      F   Pr(>F)   
Treatment  1   0.2239 0.03688 0.8425 0.005994 **
Residual  22   5.8474 0.96312                   
Total     23   6.0713 1.00000                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

attr(,"class")
[1] "pwadstrata" "list"      

Time point 1 (pre-treatment): Significant effect of Site (r2 of 0.73). Treatment is also a significant predictor (r2 of 0.04). This is something to note because this is supposed to be before treatment was applied. When permutations are restricted to within site, treatment and site remain significant. For weighted unifrac, treatment is insignificant (r2 of 0.025) whether site is controlled for or not.

#additional analysis of treatment effects when controling for site.
adonis2(dist ~  Treatment , data= sd_adonis, strata = sd_adonis$Site, permutations = 1000)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  strata 
Permutation: free
Number of permutations: 1000

adonis2(formula = dist ~ Treatment, data = sd_adonis, permutations = 1000, strata = sd_adonis$Site)
          Df SumOfSqs      R2      F  Pr(>F)  
Treatment  2   0.4285 0.05346 0.9318 0.02797 *
Residual  33   7.5869 0.94654                 
Total     35   8.0154 1.00000                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#and pairwise
pairwise.adonis2(dist ~ Treatment, strata = "Site",  data= sd_adonis, nperm = 1000)

$parent_call
[1] "dist ~ Treatment , strata = Site , permutations 1000"

$CF_vs_CON
          Df SumOfSqs     R2      F  Pr(>F)  
Treatment  1   0.2950 0.0568 1.3248 0.01199 *
Residual  22   4.8983 0.9432                 
Total     23   5.1933 1.0000                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$CF_vs_ORG
          Df SumOfSqs      R2      F  Pr(>F)  
Treatment  1   0.2518 0.05346 1.2424 0.01998 *
Residual  22   4.4592 0.94654                 
Total     23   4.7110 1.00000                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

$CON_vs_ORG
          Df SumOfSqs      R2      F Pr(>F)
Treatment  1   0.0959 0.01622 0.3628 0.7093
Residual  22   5.8163 0.98378              
Total     23   5.9122 1.00000              

attr(,"class")
[1] "pwadstrata" "list"      

Time point 2 (post-treatment): Site was again a significant predictor (r2 of 0.60) and treatment was borderline significant (r2 of 0.03). However, when the effect of site was controlled for by using strata = site, treatment becomes significant (r2 of 0.039). So when controlling for site differences, the treatments contribute to differences in microbiome structure. This was not the case for weighted unifrac dissimilarity. Treatment was not a significant predictor of community dissimilarity regardless of whether or not site differences were controlled for.

Take away: Site remains the most useful predictor of community dissimilarity. Treatment effects are only observed for non-phylogentic metrics of dissimilarity when the effect of site is controlled for.

#Visualize ordinations using two methods. CAP and NMDS

#extract data from phyloseq object and run anova.cca and drop1(cca) using vegan to examine significance of constraining variables
sd_adonis<-data.frame(sample_data(bb16S_bac_ML))
sd_adonis$Treatment <-as.factor(sd_adonis$Treatment)
sd_adonis$Site <-as.factor(sd_adonis$Site)
tab_adonis<-data.frame(otu_table(bb16S_bac_ML))
#run model and check global significance using anova.cca.
cca_mod<-cca(tab_adonis ~ Site + Timepoint + Treatment, sd_adonis)
anova.cca(cca_mod)

Permutation test for cca under reduced model
Permutation: free
Number of permutations: 999

Model: cca(formula = tab_adonis ~ Site + Timepoint + Treatment, data = sd_adonis)
         Df ChiSquare      F Pr(>F)    
Model     8    1.6817 3.0164  0.001 ***
Residual 63    4.3906                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#significant global model so lets look at individual parameters
drop1(cca_mod, test = "perm", permutations = 100)

          Df    AIC      F   Pr(>F)   
<none>       124.52                   
Site       5 133.42 3.7809 0.009901 **
Timepoint  1 126.09 3.2013 0.009901 **
Treatment  2 122.80 1.0125 0.386139   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

0#Site and time point are significant but treatment is not. We will not include treatment in our constrained ordination.

[1] 0

Constrained ordination using CAP

#Bray-Curtis
ord_cap<-ordinate(physeq = bb16S_bac_ML, method = "CAP", distance = "bray", formula = ~  Site + Timepoint )
plot_ordination(physeq = bb16S_bac_ML, ordination = ord_cap, type = "Samples", color = "Site", shape = "Timepoint" ) + theme_classic() + ggtitle("Bray-Curtis CAP")

#strongest grouping by Site. Time point also has a smaller effect and the groupings are not nearly as clear. 

ord_cap<-ordinate(physeq = bb16S_bac_ML, method = "CAP", distance = "wunifrac", formula = ~  Site + Timepoint )
plot_ordination(physeq = bb16S_bac_ML, ordination = ord_cap, type = "Samples", color = "Site", shape = "Timepoint" ) + theme_classic()  + ggtitle("Weighted Unifrac CAP")

#strongest grouping by Site. Time point also has a smaller effect and the groupings are not nearly as clear. There are several outlier samples that ordinate far to the right on the CAP1 axis. 

Unconstrained ordination using NMDS

#Bray-Curtis
ord_nmds<-ordinate(physeq = bb16S_bac_ML, method = "NMDS", distance = "bray", k = 3, trymax = 1000)

Run 0 stress 0.103921 
Run 1 stress 0.1056303 
Run 2 stress 0.1053541 
Run 3 stress 0.1150856 
Run 4 stress 0.1037439 
... New best solution
... Procrustes: rmse 0.01088478  max resid 0.08656313 
Run 5 stress 0.1059481 
Run 6 stress 0.1156754 
Run 7 stress 0.1039219 
... Procrustes: rmse 0.0108877  max resid 0.08661856 
Run 8 stress 0.1037624 
... Procrustes: rmse 0.001371454  max resid 0.01002291 
Run 9 stress 0.117372 
Run 10 stress 0.1136756 
Run 11 stress 0.1140731 
Run 12 stress 0.1037441 
... Procrustes: rmse 0.0001591934  max resid 0.0004472688 
... Similar to previous best
Run 13 stress 0.1121667 
Run 14 stress 0.1068808 
Run 15 stress 0.103921 
... Procrustes: rmse 0.01086929  max resid 0.0865921 
Run 16 stress 0.1133165 
Run 17 stress 0.1158687 
Run 18 stress 0.1056302 
Run 19 stress 0.105631 
Run 20 stress 0.1039207 
... Procrustes: rmse 0.01084584  max resid 0.08654033 
*** Solution reached

ord_nmds$stress

[1] 0.1037439

plot_ordination(physeq = bb16S_bac_ML, ordination = ord_nmds, type = "Samples", color = "Site", shape = "Timepoint" ) + theme_classic() + ggtitle("Bray-Curtis NMDS")

#strongest grouping by Site. Time point also has a smaller effect, but these groupings are more clear than in the CAP ordination. 

#weighted Unifrac
ord_nmds<-ordinate(physeq = bb16S_bac_ML, method = "NMDS", distance = "wunifrac")

Run 0 stress 0.04895024 
Run 1 stress 0.04895017 
... New best solution
... Procrustes: rmse 2.032119e-05  max resid 0.0001337503 
... Similar to previous best
Run 2 stress 0.0489501 
... New best solution
... Procrustes: rmse 0.0001782699  max resid 0.001032329 
... Similar to previous best
Run 3 stress 0.04895017 
... Procrustes: rmse 2.374303e-05  max resid 0.0001527707 
... Similar to previous best
Run 4 stress 0.04895002 
... New best solution
... Procrustes: rmse 5.859808e-05  max resid 0.0003491409 
... Similar to previous best
Run 5 stress 0.04895014 
... Procrustes: rmse 7.271456e-05  max resid 0.0004277203 
... Similar to previous best
Run 6 stress 0.04895011 
... Procrustes: rmse 9.595678e-05  max resid 0.0005658183 
... Similar to previous best
Run 7 stress 0.04895011 
... Procrustes: rmse 6.251135e-05  max resid 0.0003675656 
... Similar to previous best
Run 8 stress 0.04895005 
... Procrustes: rmse 2.790967e-05  max resid 0.0001631977 
... Similar to previous best
Run 9 stress 0.04895016 
... Procrustes: rmse 8.040287e-05  max resid 0.0004696986 
... Similar to previous best
Run 10 stress 0.04895028 
... Procrustes: rmse 0.0001411818  max resid 0.0008290002 
... Similar to previous best
Run 11 stress 0.04895012 
... Procrustes: rmse 6.358153e-05  max resid 0.0003773398 
... Similar to previous best
Run 12 stress 0.04895024 
... Procrustes: rmse 0.0001036116  max resid 0.0005910888 
... Similar to previous best
Run 13 stress 0.04895025 
... Procrustes: rmse 9.601491e-05  max resid 0.0005456503 
... Similar to previous best
Run 14 stress 0.0489502 
... Procrustes: rmse 9.320514e-05  max resid 0.000557752 
... Similar to previous best
Run 15 stress 0.04895018 
... Procrustes: rmse 8.698134e-05  max resid 0.0005129484 
... Similar to previous best
Run 16 stress 0.04895017 
... Procrustes: rmse 8.289461e-05  max resid 0.0004970631 
... Similar to previous best
Run 17 stress 0.04895011 
... Procrustes: rmse 9.887275e-05  max resid 0.0005832849 
... Similar to previous best
Run 18 stress 0.04895011 
... Procrustes: rmse 5.956548e-05  max resid 0.0003440157 
... Similar to previous best
Run 19 stress 0.04895023 
... Procrustes: rmse 0.000138979  max resid 0.0008135799 
... Similar to previous best
Run 20 stress 0.04895025 
... Procrustes: rmse 0.0001047004  max resid 0.000616414 
... Similar to previous best
*** Solution reached

ord_nmds$stress

[1] 0.04895002

plot_ordination(physeq = bb16S_bac_ML, ordination = ord_nmds, type = "Samples", color = "Site", shape = "Timepoint" ) + theme_classic() + ggtitle("Weighted Unifrac NMDS")

#strongest grouping by Site. Time point also has a smaller effect, but these groupings are more clear than in the CAP ordination. There are several outlier samples that ordinate far to the right on the NMDS1 axis. 

#Indicator taxa using Baruta Algorith Specify function for Baruta that accepts a phyloseq object.

baruta_phyloseq<-function(physeq = "Phyloseq object", GroupingVar = "Grouping variable"){
  #prep data
  df<-data.frame((otu_table(physeq)))
  md<-data.frame(sample_data(physeq))
  tax<-data.frame(tax_table(physeq))
  rownames(tax)<-str_replace(rownames(tax), pattern = "[[.]]", replacement = "=")
  #run simper
 
                  boruta_layer<-Boruta(df , factor(md[,GroupingVar]), 
                      doTrace = 1, ntree = 1000, maxRuns = 1000) # adjust factor!!!!

                  B1<-data.frame(boruta_layer$finalDecision)  
                  B2<-data.frame(apply(boruta_layer$ImpHistory, 2, mean))
                  names(B1)<-"Layer_decision"
                  names(B2)<-"Layer_importance"
                  BB1<-merge(B1, B2, by="row.names", all.x=TRUE)
                  #sorting the results
                  BB1_sort<-BB1[order(BB1$Layer_importance, decreasing = TRUE),] ## this is to give an overview and choose cutoff
                  #picking up important with cutoff
                  BB1_sort_cutoff<-filter(BB1, Layer_importance >= 2) ## number is your importance cutoff

                        #picking up OTUs from full otutable
                        signOTU<-as.matrix(BB1_sort_cutoff$Row.names) #vector with OTUs
                        signOTU<-str_replace(signOTU, pattern = "[[.]]", replacement = ".")
                        #my_taxa
                        tax$OTU <-rownames(tax)
                        signOTU<-str_replace(signOTU, pattern = "[[.]]", replacement = "=")
                        signmy_taxa<-filter(tax, OTU %in% signOTU) 
                        signmy_taxa<-as.data.frame(signmy_taxa)
                        taxa_keep<-signmy_taxa$OTU
return(taxa_keep)
}

Run Baruta algorith with full data set. We are looking for taxa that explain differences among sites.

baruta_Site_indicators<-baruta_phyloseq(physeq = bb16S_bac_ML, GroupingVar = "Site")
baruta_Site_indicators_ps_object<-subset_taxa(bb16S_bac_ML, taxa_names(bb16S_bac_ML) %in% baruta_Site_indicators)

https://jokergoo.github.io/ComplexHeatmap-reference/book/ Heatmap to visualize differnces with all samples together.

#create subset phyloseq object and pull out taxonomic info
baruta_Site_indicators_ps_object_tax<-tax_table(baruta_Site_indicators_ps_object)
baruta_Site_indicators_ps_object_tax<-data.frame(baruta_Site_indicators_ps_object_tax)
#rename NA as unknown
baruta_Site_indicators_ps_object_tax$Class[is.na(baruta_Site_indicators_ps_object_tax$Class)]<-"Unknown"
baruta_Site_indicators_ps_object_tax$Family[is.na(baruta_Site_indicators_ps_object_tax$Family)]<-"Unknown"
baruta_Site_indicators_ps_object@tax_table<-tax_table(baruta_Site_indicators_ps_object_tax)
fam_names<-data.frame(gsub("column", "", tax_table(baruta_Site_indicators_ps_object)[,3]))

#create color pallete. This goes from green to white depending on z.score of individual taxa
col_fun<-colorRamp2(c(-2, 0, 2), c("green", "white", "red"))
#Heatmap with taxa names
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = tax_table(baruta_Site_indicators_ps_object)[,4], row_dend_reorder =TRUE, col= col_fun)
#heatmap with otuIDs
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = taxa_names(baruta_Site_indicators_ps_object), row_dend_reorder =TRUE, col= col_fun)

Based off of initial figure, it looks like there is a strong grouping based on Baruta ID’ed taxa, and there is pattern for every other row. This indicates differences in time points within each sites. We might want to split by time point and then repeat analysis.

Since we observe what appears to be a time point signature, we can split by time point and rerun the Baruta algorithm.

#pre-treatment samples
bb16S_bac_ML_A<-subset_samples(bb16S_bac_ML, Timepoint == "A")
#run algo and subset phyloseq object
baruta_Site_indicators_A<-baruta_phyloseq(physeq = bb16S_bac_ML_A, GroupingVar = "Site")
baruta_Site_indicators_ps_object_A<-subset_taxa(bb16S_bac_ML_A, taxa_names(bb16S_bac_ML_A) %in% baruta_Site_indicators_A)
# pull out taxonomic info
baruta_Site_indicators_ps_object_tax_A<-tax_table(baruta_Site_indicators_ps_object_A)
baruta_Site_indicators_ps_object_tax_A<-data.frame(baruta_Site_indicators_ps_object_tax_A)
#rename NA as unknown
baruta_Site_indicators_ps_object_tax_A$Class[is.na(baruta_Site_indicators_ps_object_tax_A$Class)]<-"Unknown"
baruta_Site_indicators_ps_object_tax_A$Family[is.na(baruta_Site_indicators_ps_object_tax_A$Family)]<-"Unknown"
#pull out taxonomy table for naming purposes
baruta_Site_indicators_ps_object_A@tax_table<-tax_table(baruta_Site_indicators_ps_object_tax_A)
fam_names<-data.frame(gsub("column", "", tax_table(baruta_Site_indicators_ps_object_A)[,3]))

#heatmap with taxa names
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_A, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = tax_table(baruta_Site_indicators_ps_object_A)[,4], row_dend_reorder =TRUE, col= col_fun)
#heatmap with otuIDs
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_A, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = taxa_names(baruta_Site_indicators_ps_object_A), row_dend_reorder =TRUE, col= col_fun)

Strong grouping of sites based on Baruta ID’ed taxa at the first sampling.

#post-treatment samples
bb16S_bac_ML_B<-subset_samples(bb16S_bac_ML, Timepoint == "B")
#run algo and subset phyloseq object
baruta_Site_indicators_B<-baruta_phyloseq(physeq = bb16S_bac_ML_B, GroupingVar = "Site")
baruta_Site_indicators_ps_object_B<-subset_taxa(bb16S_bac_ML_B, taxa_names(bb16S_bac_ML_B) %in% baruta_Site_indicators_B)
#extract taxonomy info and rename NAs as unknown
baruta_Site_indicators_ps_object_tax_B<-tax_table(baruta_Site_indicators_ps_object_B)
baruta_Site_indicators_ps_object_tax_B<-data.frame(baruta_Site_indicators_ps_object_tax_B)
baruta_Site_indicators_ps_object_tax_B$Class[is.na(baruta_Site_indicators_ps_object_tax_B$Class)]<-"Unknown"
baruta_Site_indicators_ps_object_tax_B$Family[is.na(baruta_Site_indicators_ps_object_tax_B$Family)]<-"Unknown"
baruta_Site_indicators_ps_object_B@tax_table<-tax_table(baruta_Site_indicators_ps_object_tax_B)
fam_names<-data.frame(gsub("column", "", tax_table(baruta_Site_indicators_ps_object_B)[,3]))

#heatmap with taxa names
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_B, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = tax_table(baruta_Site_indicators_ps_object_B)[,4], row_dend_reorder =TRUE, col= col_fun)
#heatmap with OTUids
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_B, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = taxa_names(baruta_Site_indicators_ps_object_B), row_dend_reorder =TRUE, col= col_fun)

Again, strong groupings of sites by Baruta ID’ed taxa at time point 2.

Next, we will look for taxa that consistently explain the differences among treatments following application. Again we will use the Baruta algorithm, but this time, we will specify the grouping variable as Treatment.

Merge by treatments to see if there are consistencies across treatments. No significant indicators.

#subset to post-treatment samples only
bb16S_bac_ML_B<-subset_samples(bb16S_bac_ML, Timepoint == "B")
#summarize by Treatment_Site
bb16S_bac_ML_B_merged<-merge_samples(bb16S_bac_ML_B, "Treatment")
sample_data(bb16S_bac_ML_B_merged)$Treatment<-rownames(sample_data(bb16S_bac_ML_B_merged))
#run algo and filter phyloseq object
bb16S_bac_ML_B_merged<-prune_taxa(taxa_sums(bb16S_bac_ML_B_merged) > 0, bb16S_bac_ML_B_merged)
baruta_Site_indicators_B_merged<-baruta_phyloseq(physeq = bb16S_bac_ML_B_merged, GroupingVar = "Treatment")
#no taxa indicative of treatment x site interaction.
#Stop this inquiry here. Do not run below because no taxa are identified. 

baruta_Site_indicators_ps_object_B_merged<-subset_taxa(bb16S_bac_ML_B_merged, taxa_names(bb16S_bac_ML_B_merged) %in% baruta_Site_indicators_B_merged)
#extract taxonomy and rename NA to unknown. 
baruta_Site_indicators_ps_object_tax_B_merged<-tax_table(baruta_Site_indicators_ps_object_B_merged)
baruta_Site_indicators_ps_object_tax_B_merged<-data.frame(baruta_Site_indicators_ps_object_tax_B_merged)
baruta_Site_indicators_ps_object_tax_B_merged$Class[is.na(baruta_Site_indicators_ps_object_tax_B_merged$Class)]<-"Unknown"
baruta_Site_indicators_ps_object_tax_B_merged$Family[is.na(baruta_Site_indicators_ps_object_tax_B_merged$Family)]<-"Unknown"
baruta_Site_indicators_ps_object_B_merged@tax_table<-tax_table(baruta_Site_indicators_ps_object_tax_B_merged)
fam_names<-data.frame(gsub("column", "", tax_table(baruta_Site_indicators_ps_object_B_merged)[,3]))
#heatmap with taxa names
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_B_merged, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = tax_table(baruta_Site_indicators_ps_object_B)[,4], row_dend_reorder =TRUE, col= col_fun)
#heatmap with OTUids
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_B_merged, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = taxa_names(baruta_Site_indicators_ps_object_B), row_dend_reorder =TRUE, col= col_fun)

Merge by treatment site and rerun algo. No significant indicators.

#subset to post-treatment samples only
bb16S_bac_ML_B<-subset_samples(bb16S_bac_ML, Timepoint == "B")
#summarize by Treatment_Site
bb16S_bac_ML_B_merged<-merge_samples(bb16S_bac_ML_B, "Treatment_Site")
sample_data(bb16S_bac_ML_B_merged)$Treatment_Site<-rownames(sample_data(bb16S_bac_ML_B_merged))
#run algo and filter phyloseq object
bb16S_bac_ML_B_merged<-prune_taxa(taxa_sums(bb16S_bac_ML_B_merged) > 0, bb16S_bac_ML_B_merged)
baruta_Site_indicators_B_merged<-baruta_phyloseq(physeq = bb16S_bac_ML_B_merged, GroupingVar = "Treatment_Site")
#no taxa indicative of treatment x site interaction.
#Stop this inquiry here. Do not run below becasue no taxa are identified. 
baruta_Site_indicators_ps_object_B_merged<-subset_taxa(bb16S_bac_ML_B_merged, taxa_names(bb16S_bac_ML_B_merged) %in% baruta_Site_indicators_B_merged)
#extract taxonomy and rename NA to unknown. 
baruta_Site_indicators_ps_object_tax_B_merged<-tax_table(baruta_Site_indicators_ps_object_B_merged)
baruta_Site_indicators_ps_object_tax_B_merged<-data.frame(baruta_Site_indicators_ps_object_tax_B_merged)
baruta_Site_indicators_ps_object_tax_B_merged$Class[is.na(baruta_Site_indicators_ps_object_tax_B_merged$Class)]<-"Unknown"
baruta_Site_indicators_ps_object_tax_B_merged$Family[is.na(baruta_Site_indicators_ps_object_tax_B_merged$Family)]<-"Unknown"
baruta_Site_indicators_ps_object_B_merged@tax_table<-tax_table(baruta_Site_indicators_ps_object_tax_B_merged)
fam_names<-data.frame(gsub("column", "", tax_table(baruta_Site_indicators_ps_object_B_merged)[,3]))
#heatmap with taxa names
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_B_merged, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = tax_table(baruta_Site_indicators_ps_object_B)[,4], row_dend_reorder =TRUE, col= col_fun)
#heatmap with OTUids
micrUBIfuns::plot_taxa_heatmap(baruta_Site_indicators_ps_object_B_merged, rm_na = F,scale_by = "taxa", cluster_rows=F, cluster_columns = T ,column_labels = taxa_names(baruta_Site_indicators_ps_object_B), row_dend_reorder =TRUE, col= col_fun)

Exploratory tree of indicators from initial Baruta and then split by timpeoints.

plot_tree(baruta_Site_indicators_ps_object, color="Site", shape = "Timepoint", label.tips = "taxa_names")

Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'physeq' in selecting a method for function 'phy_tree': object 'baruta_Site_indicators_ps_object' not found

#create super long color vector
col_vector <- c("#000000", "#FFFF00", "#1CE6FF", "#FF34FF", "#FF4A46", "#008941", "#006FA6", "#A30059",
        "#FFDBE5", "#7A4900", "#0000A6", "#63FFAC", "#B79762", "#004D43", "#8FB0FF", "#997D87",
        "#5A0007", "#809693", "#FEFFE6", "#1B4400", "#4FC601", "#3B5DFF", "#4A3B53", "#FF2F80",
        "#61615A", "#BA0900", "#6B7900", "#00C2A0", "#FFAA92", "#FF90C9", "#B903AA", "#D16100",
        "#DDEFFF", "#000035", "#7B4F4B", "#A1C299", "#300018", "#0AA6D8", "#013349", "#00846F",
        "#372101", "#FFB500", "#C2FFED", "#A079BF", "#CC0744", "#C0B9B2", "#C2FF99", "#001E09",
        "#00489C", "#6F0062", "#0CBD66", "#EEC3FF", "#456D75", "#B77B68", "#7A87A1", "#788D66",
        "#885578", "#FAD09F", "#FF8A9A", "#D157A0", "#BEC459", "#456648", "#0086ED", "#886F4C",
        
        "#34362D", "#B4A8BD", "#00A6AA", "#452C2C", "#636375", "#A3C8C9", "#FF913F", "#938A81",
        "#575329", "#00FECF", "#B05B6F", "#8CD0FF", "#3B9700", "#04F757", "#C8A1A1", "#1E6E00",
        "#7900D7", "#A77500", "#6367A9", "#A05837", "#6B002C", "#772600", "#D790FF", "#9B9700",
        "#549E79", "#FFF69F", "#201625", "#72418F", "#BC23FF", "#99ADC0", "#3A2465", "#922329",
        "#5B4534", "#FDE8DC", "#404E55", "#0089A3", "#CB7E98", "#A4E804", "#324E72", "#6A3A4C",
        "#83AB58", "#001C1E", "#D1F7CE", "#004B28", "#C8D0F6", "#A3A489", "#806C66", "#222800",
        "#BF5650", "#E83000", "#66796D", "#DA007C", "#FF1A59", "#8ADBB4", "#1E0200", "#5B4E51",
        "#C895C5", "#320033", "#FF6832", "#66E1D3", "#CFCDAC", "#D0AC94", "#7ED379", "#012C58",
        
        "#7A7BFF", "#D68E01", "#353339", "#78AFA1", "#FEB2C6", "#75797C", "#837393", "#943A4D",
        "#B5F4FF", "#D2DCD5", "#9556BD", "#6A714A", "#001325", "#02525F", "#0AA3F7", "#E98176",
        "#DBD5DD", "#5EBCD1", "#3D4F44", "#7E6405", "#02684E", "#962B75", "#8D8546", "#9695C5",
        "#E773CE", "#D86A78", "#3E89BE", "#CA834E", "#518A87", "#5B113C", "#55813B", "#E704C4",
        "#00005F", "#A97399", "#4B8160", "#59738A", "#FF5DA7", "#F7C9BF", "#643127", "#513A01",
        "#6B94AA", "#51A058", "#A45B02", "#1D1702", "#E20027", "#E7AB63", "#4C6001", "#9C6966",
        "#64547B", "#97979E", "#006A66", "#391406", "#F4D749", "#0045D2", "#006C31", "#DDB6D0",
        "#7C6571", "#9FB2A4", "#00D891", "#15A08A", "#BC65E9", "#FFFFFE", "#C6DC99", "#203B3C",
        "#671190", "#6B3A64", "#F5E1FF", "#FFA0F2", "#CCAA35", "#374527", "#8BB400", "#797868",
        "#C6005A", "#3B000A", "#C86240", "#29607C", "#402334", "#7D5A44", "#CCB87C", "#B88183",
        "#AA5199", "#B5D6C3", "#A38469", "#9F94F0", "#A74571", "#B894A6", "#71BB8C", "#00B433",
        "#789EC9", "#6D80BA", "#953F00", "#5EFF03", "#E4FFFC", "#1BE177", "#BCB1E5", "#76912F",
        "#003109", "#0060CD", "#D20096", "#895563", "#29201D", "#5B3213", "#A76F42", "#89412E",
        "#1A3A2A", "#494B5A", "#A88C85", "#F4ABAA", "#A3F3AB", "#00C6C8", "#EA8B66", "#958A9F",
        "#BDC9D2", "#9FA064", "#BE4700", "#658188", "#83A485", "#453C23", "#47675D", "#3A3F00",
        "#061203", "#DFFB71", "#868E7E", "#98D058", "#6C8F7D", "#D7BFC2", "#3C3E6E", "#D83D66",
        "#2F5D9B", "#6C5E46", "#D25B88", "#5B656C", "#00B57F", "#545C46", "#866097", "#365D25",
        "#252F99", "#00CCFF", "#674E60", "#FC009C", "#92896B")

#Core and unique taxa

#core adonis
tab_adonis<-data.frame(otu_table(core_asv))
sd_adonis<-data.frame(sample_data(core_asv))
sd_adonis$Treatment <-as.factor(sd_adonis$Treatment)
sd_adonis$Site <-as.factor(sd_adonis$Site)
dist<-phyloseq::distance(core_asv, method =  "bray")

#full model with no strata argument. Essentially the effect of each main effect without controlling for any potential groups. 
adonis2(dist ~  Site * Treatment * Timepoint, data= sd_adonis, permutations = 1000)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 1000

adonis2(formula = dist ~ Site * Treatment * Timepoint, data = sd_adonis, permutations = 1000)
                         Df SumOfSqs      R2       F   Pr(>F)    
Site                      5   8.7827 0.55648 26.0777 0.000999 ***
Treatment                 2   0.4337 0.02748  3.2192 0.009990 ** 
Timepoint                 1   0.7005 0.04438 10.3998 0.000999 ***
Site:Treatment           10   1.0815 0.06853  1.6056 0.044955 *  
Site:Timepoint            5   1.3256 0.08399  3.9361 0.000999 ***
Treatment:Timepoint       2   0.2154 0.01365  1.5990 0.163836    
Site:Treatment:Timepoint 10   0.8183 0.05185  1.2149 0.219780    
Residual                 36   2.4249 0.15364                     
Total                    71  15.7827 1.00000                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Core take aways: The core was defined as those taxa occuring in atleast 90% of samples and having a relative abudnance of greater than 0.00001. When collapsed at the genus level, there are 33 members in the core. Using ASVs, the core consists of 49 members. Looking at the phylogenetic trees, the members of Rickettsiales condense into a single tip at the genus level. Blast-ing this sequence returns plant DNA. I wonder if this is a common plant visited by these bees? Either way, the core accoutns for a mean 27% +/- 0.004 (sd) of the total reads in each sample.

Taxonomic composition of core: Phylum: Actinobacteriota Bacteroidota Firmicutes Proteobacteria 2 2 8 19

Class: Actinobacteria Alphaproteobacteria Bacilli Bacteroidia Clostridia Gammaproteobacteria 2 8 6 2 2 11

non-core taxa

tab_adonis<-data.frame(otu_table(non_core_asv))
sd_adonis<-data.frame(sample_data(non_core_asv))
sd_adonis$Treatment <-as.factor(sd_adonis$Treatment)
sd_adonis$Site <-as.factor(sd_adonis$Site)
dist<-phyloseq::distance(non_core_asv, method =  "bray")

#full model with no strata argument. Essentially the effect of each main effect without controlling for any potential groups. 
adonis2(dist ~  Site * Treatment * Timepoint, data= sd_adonis, permutations = 1000)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 1000

adonis2(formula = dist ~ Site * Treatment * Timepoint, data = sd_adonis, permutations = 1000)
                         Df SumOfSqs      R2      F   Pr(>F)    
Site                      5   6.2949 0.20248 4.5047 0.000999 ***
Treatment                 2   0.7228 0.02325 1.2932 0.042957 *  
Timepoint                 1   2.1193 0.06817 7.5832 0.000999 ***
Site:Treatment           10   3.3110 0.10650 1.1847 0.018981 *  
Site:Timepoint            5   4.8300 0.15536 3.4564 0.000999 ***
Treatment:Timepoint       2   0.6569 0.02113 1.1751 0.133866    
Site:Treatment:Timepoint 10   3.0934 0.09950 1.1069 0.106893    
Residual                 36  10.0613 0.32362                    
Total                    71  31.0897 1.00000                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Determine those taxa unique to a site at Genus level
unique_genera_to_region_tax<-unique_taxa(phyloseq_obj = bb16S_bac_rarefy_tax_genus, treatment = "Region")
unique_genera_to_site_tax<-unique_taxa(phyloseq_obj = bb16S_bac_rarefy_tax_genus, treatment = "Site")
#asv level
unique_asv_to_region_tax<-unique_taxa(phyloseq_obj = bb16S_bac_rarefy, treatment = "Region")
unique_asv_to_site_tax<-unique_taxa(phyloseq_obj = bb16S_bac_rarefy, treatment = "Site")

#genera
#186 genera are unique to one region. 176 genera are unique to a site. This means that only 10 Taxa unique to a region are not shared for the sites within that region. 
unique_genera_to_site_tax_ps<-subset_taxa(bb16S_bac_rarefy, taxa_names(bb16S_bac_rarefy) %in% unlist(unique_genera_to_site_tax))
unique_genera_to_region_tax_ps<-subset_taxa(bb16S_bac_rarefy, taxa_names(bb16S_bac_rarefy) %in% unlist(unique_genera_to_region_tax))

#asv
#1510 asvs are unique to one region. 1484 asvs are unique to a site. This means that only 26 asvs unique to a region are not shared for the sites within that region. 
unique_asv_to_site_tax_ps<-subset_taxa(bb16S_bac_rarefy, taxa_names(bb16S_bac_rarefy) %in% unlist(unique_asv_to_site_tax))
unique_asv_to_region_tax_ps<-subset_taxa(bb16S_bac_rarefy, taxa_names(bb16S_bac_rarefy) %in% unlist(unique_asv_to_region_tax))

plot_bar(unique_genera_to_site_tax_ps, x = "Region_Site", fill = "Genus") + scale_fill_manual(values = col_vector) +geom_bar(stat="identity")
plot_bar(unique_genera_to_region_tax_ps, x = "Region", fill = "Genus") + scale_fill_manual(values = col_vector) +geom_bar(stat="identity")

#plot treees
plot_tree(unique_genera_to_site_tax_ps, color = "Site", shape = "Timepoint", label.tips = "Order", text.size = 3 )
plot_tree(unique_genera_to_region_tax_ps, color = "Site", shape = "Timepoint", label.tips = "Order", text.size = 3 )


#percentage of reads accounted for by core. 
summary(sample_sums(unique_genera_to_site_tax_ps)/25000)
summary(sample_sums(unique_genera_to_region_tax_ps)/25000)


#taxonomic composition
table(tax_table(core_genus)[,2])
table(tax_table(core_genus)[,3])


#pie charts of endemic taxa taxonomy
#Filter to regions and create pie charts for each. I've chosen to plot regions becasue its a bit simplier but we could quickly change this to sites. 
#CP
CP<-subset_samples(unique_genera_to_site_tax_ps, Region == "CP")
CP<-filter_taxa(CP, function(x) sum(x) >= 1 , TRUE)
pie_data_CP<-data.frame(table(tax_table(CP)[,2]))
names(pie_data_CP) <- c("Class", "Frequency")
ggplot(pie_data_CP, aes(x="", y=Frequency, fill=Class)) +
  geom_bar(stat="identity", width=1, color="white") +
  coord_polar("y", start=0) +
  theme_void()   + scale_fill_manual(values = col_vector) + ggtitle("CP endemic genera at Phylum")

#NE
NE<-subset_samples(unique_genera_to_site_tax_ps, Region == "NE")
NE<-filter_taxa(NE, function(x) sum(x) >= 1 , TRUE)
pie_data_NE<-data.frame(table(tax_table(NE)[,2]))
names(pie_data_NE) <- c("Class", "Frequency")
ggplot(pie_data_NE, aes(x="", y=Frequency, fill=Class)) +
  geom_bar(stat="identity", width=1, color="white") +
  coord_polar("y", start=0) +
  theme_void()   + scale_fill_manual(values = col_vector) + ggtitle("NE endemic genera at Phylum")

#SE
SE<-subset_samples(unique_genera_to_site_tax_ps, Region == "SE")
SE<-filter_taxa(SE, function(x) sum(x) >= 1 , TRUE)
pie_data_SE<-data.frame(table(tax_table(SE)[,2]))
names(pie_data_SE) <- c("Class", "Frequency")
ggplot(pie_data_SE, aes(x="", y=Frequency, fill=Class)) +
  geom_bar(stat="identity", width=1, color="white") +
  coord_polar("y", start=0) +
  theme_void()   + scale_fill_manual(values = col_vector) + ggtitle("SE endemic genera at Phylum")

#WV
WV<-subset_samples(unique_genera_to_site_tax_ps, Region == "WV")
WV<-filter_taxa(WV, function(x) sum(x) >= 1 , TRUE)
pie_data_WV<-data.frame(table(tax_table(WV)[,2]))
names(pie_data_WV) <- c("Class", "Frequency")
ggplot(pie_data_WV, aes(x="", y=Frequency, fill=Class)) +
  geom_bar(stat="identity", width=1, color="white") +
  coord_polar("y", start=0) +
  theme_void()   + scale_fill_manual(values = col_vector) + ggtitle("WV endemic genera at Phylum")

Endemic taxa take aways: 186 genera are unique to one region. 176 genera are unique to a site. This means that only 10 Taxa unique to a region are not shared for the sites within that region. Endemic taxa are currently defined as those occurring only at a single site. This creates its own set of issues in that many taxa are only found in a single sample. I’ve reached out to the PhyloSmith package creator to ask about augmenting this funciton to relax this a bit to include those taxa occurring in two or three sites. Endemic taxa are much less abundant and account for a max of 0.0086 % of the reads at a site (mean 0.00064%).

Alternate approach to core and endemic taxa based on phylognetic groupings as opposed to collpasing by taxonomy. Shared and unique taxa among sites by phylogeny

#cut tree at two differnt levels for exploratory purposes.
bb16S_bac_rarefy_tip_merged0_1<-tip_glom(physeq = bb16S_bac_rarefy, h = 0.1)
plot_tree(bb16S_bac_rarefy_tip_merged0_1, label.tips = "Family", text.size = 3, color = "Region")
#275 members
bb16S_bac_rarefy_tip_merged0_2<-tip_glom(physeq = bb16S_bac_rarefy, h = 0.2)
plot_tree(bb16S_bac_rarefy_tip_merged0_2, label.tips = "Family", text.size = 3, color = "Region")
#core - 24 members in core 
core_genus_tip<-taxa_core(phyloseq_obj = bb16S_bac_rarefy_tip_merged0_2, treatment = "Site", frequency = 0.90, abundance_threshold = 0.00001 ) 
#plot tip core
plot_tree(core_genus_tip, color = "Site", shape = "Timepoint", label.tips = "Family", text.size = 3 )


#Determine those taxa unique to a site in 0.2 cut tree
unique_genera_to_region_tip<-unique_taxa(phyloseq_obj = bb16S_bac_rarefy_tip_merged0_2, treatment = "Region")
unique_genera_to_site_tip<-unique_taxa(phyloseq_obj = bb16S_bac_rarefy_tip_merged0_2, treatment = "Site")

# 137 unique to a site and 144 unique to a region. 
unique_genera_to_site_tip_ps<-subset_taxa(bb16S_bac_rarefy, taxa_names(bb16S_bac_rarefy) %in% unlist(unique_genera_to_site_tip))
plot_tree(unique_genera_to_site_tip_ps, color = "Site", shape = "Timepoint", label.tips = "Order", text.size = 3 )
unique_genera_to_region_tip_ps<-subset_taxa(bb16S_bac_rarefy, taxa_names(bb16S_bac_rarefy) %in% unlist(unique_genera_to_region_tip))
plot_tree(unique_genera_to_region_tip_ps, color = "Site", shape = "Timepoint", label.tips = "Order", text.size = 3 )

Relevant bee bread microbiome literature:

DOI: https://doi.org/10.1128/genomeA.00247-18